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The most common way to simplify extensive Monte-Carlo simulations of air showers is to use 
the thinning approximation. We study its effect on the physical parameters reconstructed from 
simulated showers. To this end, we have created a library of showers simulated without thinning 
with energies from 10 17 eV to 10 18 eV, various zenith angles and primaries. This library is pub- 
licly available. Physically interesting applications of the showers simulated without thinning are 
discussed. Observables reconstructed from these showers are compared to those obtained with the 
thinning approximation. The amount of artificial fluctuations introduced by thinning is estimated. 
A simple method, multisampling, is suggested which results in a controllable suppression of artificial 
fluctuations and at the same time requires less demanding computational resources as compared to 
the usual thinning. 

PACS numbers: 98.70.Sa, 96.50.sbe, 96.50.sd 



I. INTRODUCTION 

Experimental information about cosmic particles at 
very high energies is obtained through the study of at- 
mospheric showers induced by these particles and is 
hence indirect. A necessary ingredient of these studies 
is therefore a good understanding of a shower initiated 
by a primary particle with given parameters. Since the 
shower development is a complicated random process, the 
Monte-Carlo simulations are often used to model atmo- 
spheric showers [18|. Physical parameters are then recon- 
structed from the simulations and compared to real data. 

At very high energies, however, the number of parti- 
cles in a shower is so large that the simulations start to 
require unrealistic computer resources. Among several 
ways to simplify the problem and to reduce the com- 
putational time, the thinning approximation [2] is cur- 
rently the most popular one. Its key idea is to track only 
a representative set of particles; while very efficient in 
calculations and providing correct values of observables 
on average, this method introduces artificial fluctuations 
because the number of tracked particles is reduced by 
several orders of magnitude. These artificial fluctuations 
mix with natural ones and therefore reduce the precision 
of determination of physical parameters. 

The standard approach to account for natural fluctu- 
ations in the air-shower simulations is to fix all shower 
parameters and to simulate sufficient number of artifi- 
cial showers. Technically, these showers differ by initial 
random seed numbers. All interactions in a simulated 
shower are fixed by these numbers for a given thinning 
level. Random variations of these numbers result in a 
plethora of possible interaction patterns which end up in 
a distribution of an observable quantity of interest calcu- 
lated for the showers with exactly the same initial phys- 
ical parameters. This distribution thus intends to rep- 
resent intrinsic fluctuations in the shower development. 
Both the central value and the width of this distribution 



are important for physical applications. 

In practice, however, the width of the distribution 
arises from two sources: physical fluctuations and ar- 
tificial fluctuations introduced by thinning. To obtain 
the physical width alone, one should in principle perform 
simulations without thinning. This is hardly possible for 
the highest energies at the current level of computational 
techniques since one often needs to simulate thousands 
of events for a typical study. 

The aim of the present work is to estimate the relative 
size of these artificial fluctuations (for the first time it is 
done by direct comparison of showers simulated with and 
without thinning) and to develop an efficient resource- 
saving method to suppress them in realistic calculations. 

In SecHH we start with a description (Sec. Ill A)) of the 
standard thinning algorithm and explain why its use in- 
troduces additional fluctuations. Then, we briefly recall, 
in Sec. Ill B[ conventional approaches to avoid or sup- 
press these fluctuations. Sec. Ill Cl describes the library of 
showers simulated without thinning for this study. This 
library is publicly available. Sec. IIIII is devoted to a 
quantitative study of the artificial fluctuations. A new 
method, multisampling, which allows one to suppress ef- 
ficiently these unphysical fluctuations without invoking 
extensive computer resources, is suggested and discussed 
in Sec. IIVI Sec. [V] contains the discussion of the method 
and our conclusions. 



II. THINNING APPROXIMATION AND 
BEYOND 

A. Standard thinning 

The number of particles in an extended air shower 
(EAS), and hence the CPU time and disk space required 
for its full simulation, scales with the energy of the pri- 
mary particle. At energies in excess of 10 17 eV, the num- 
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ber of particles of kinetic energy above 100 MeV at the 
ground level exceeds 10 8 and the time required to simu- 
late such a shower at a computer with a few-GHz CPU 
is of order of several days. A typical vertical shower in- 
duced by a hadron of 10 18 eV requires about 100 Gb of 
disk space and a month of CPU time. Modelling indi- 
vidual showers with incident energies of about 10 20 eV is 
at the limit of realistic capabilities of modern computers; 
meanwhile one needs thousands of simulated showers for 
comparison with experimental data. 

As a result of a full simulation of a shower, one ob- 
tains the list of all particles at the ground level. This 
information is redundant for many practical purposes. 
Real ground-based experiments detect only a small frac- 
tion of these particles, so for calculating average particle 
densities one does not need to know precise coordinates 
and energies of all particles. In the thinning approxima- 
tion 0, HJ, groups of particles are replaced by effective 
representative particles with weights. 

Let us briefly recall how the thinning approximation 
works (see e.g. Ref. Q for a detailed discussion). Denote 
the primary energy by Eq and introduce a parameter e 
called the thinning level. For each subsequent interac- 
tion, consider the energies Ej of the secondary particles 
created in this interaction. If the condition 

J2Ej<eE (1) 

is satisfied, then the method prescribes to keep one of 
the secondary particles and to discard the others. The 
probability to keep the ith particle pi is proportional to 
its energy, 

p i =E i /Y l E j . 
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To the selected particle, the weight Wi — Wo/pi is as- 
signed, where wq is the weight of the initial particle of 
this interaction (wq = 1 for the particle which initiated 
the shower). 

If the condition (TT|) is not satisfied, then the so-called 
statistical thinning operates: among the secondary par- 
ticles, a subsample of ones with energies Ey < eE$ is 
considered and (one or more) effective particles are se- 
lected with probabilities 

_ Ei> 

mm(eE ,J2j' E j') 

These particles, to which the weights uv = wo/pi> are 
assigned, are kept for further simulations together with 
original particles which have had energies Ey > eE$. 

For useful values of e, the number of particles tracked 
is reduced by a factor of 10 3 - 10 6 . For a random pro- 
cess, this change in the number of particles (and con- 
sequently, in the number of interactions) results in the 
increase of fluctuations compared to the fully simulated 
process. This means that a part of fluctuations in the 



development of a shower simulated with thinning is arti- 
ficial, that is it is present neither in the full shower sim- 
ulated with e = nor in a real EAS. For a number of ap- 
plications, these fluctuations are undesirable and should 
be suppressed or at least brought under control. 



B. Standard methods to suppress fluctuations 

In the framework of the thinning method, the fluctua- 
tions are effectively suppressed by introducing the upper 
limit on the weight factor [4]. The number of "real" 
particles tracked is thus enlarged. Maximal weights for 
hadrons and for electromagnetic particles may be as- 
signed in different ways. For a given problem, the optimal 
values of the maximal weights may be found in order to 
minimise the ratio of the size of artificial fluctuations to 
the computational time. In what follows, when we refer 
to the thinning with weights limitation, we will use the 
maximal weights optimised in Ref. [3| for the calculation 
of the particle density. 

The optimal values of parameters of thinning proce- 
dure may depend on the interaction models adopted in 
simulations for a given problem. In principle, the weights 
should be optimized for each combination of the models 
(which are updated every few years) and for each partic- 
ular task (different observables, primaries, energies, etc.). 
However, this optimization requires a dedicated time con- 
suming study in each case. We suggest another approach 
to the problem in Sec. IIVI 



C. A library of showers simulated without thinning 

We have performed simulations of air showers with- 
out thinning by making use of the CORSIKA simu- 
lation code [f|. For different showers, we have used 
QGSJET 01C @ and QGSJET 11-03 [7] as high-energy 
and GHEISHA 2002d [8| as low-energy hadronic inter- 
action models. Currently, the library contains about 40 
showers induced by primary protons, gamma-rays and 
iron nuclei with energies between 10 17 eV and 10 18 eV 
and zenith angles between 0° and 45°. The show- 
ers have been simulated for the observational condi- 
tions (atmospheric depth and geomagnetic field) of ei- 
ther AGASA [9] or the Telescope Array [Io| experi- 
ments. The shower library 11] is publicly available at 
http : / /l ivni . inr . ac . ru| Detailed information about 
input parameters used for the simulation of each shower 
is available from the library website together with full 
output files. The access to the data files is provided 
freely upon request. For users not familiar with the COR- 
SIKA output format, a "Datafile reading programming 
manual" is given, containing a working example in CH — h 
Free access to the computational resources of the server 
is provided to avoid lengthy copying of the output files 
(some of which exceed 100 Gb in size). An access request 
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form along with conditions of the usage of the library are 
available from the library website. 

Given the amount of computing resources required for 
simulation, each shower simulated without thinning is 
valuable. We hope that the open library would be useful 
in studies of various physical problems, notably facing the 
improved precision of modern experiments which often 
exceeds the precision of simulations. The library is being 
continuously extended; we plan to supplement it with 
showers of higher energies in the future. 



III. SIZE OF ARTIFICIAL FLUCTUATIONS 
DUE TO THINNING 

A. Shower-by-shower comparison 

With a library of showers simulated without thinning, 
the comparison of the observables reconstructed from 
showers with and without thinning is possible. This al- 
lows one to estimate the effect of the approximation. To 
do that, for each shower without thinning (e = 0) we 
have simulated a number of showers with different thin- 
ning levels (e ^ 0). All initial parameters (including the 
random seed numbers) were kept the same as in the e = 
simulation, which enabled us to reproduce exactly the 
same first interaction in the entire set of showers. Three 
important observables — the signal density at 600 m from 
the shower axis 5(600), the muon density at 1000 m from 
the axis ^(lOOO), and the depth of the maximal shower 
development A max — were reconstructed for each of the 
showers following the data-processing operation adopted 
by the AG AS A experiment fl9j . The detector response 
was calculated with the help of GEANT simulations in 
Ref. p3. 5(600) and ^(1000) were obtained by fitting 
the correspondin g de nsity at the ground level with empir- 
ical formulae fl3l Il4j . For fitting purpose the density was 
binned into 50m-width rings centered at the shower axis. 
I max was obtained by fitting the longitudinal shower pro- 
file with the empirical Gaisser-Hillas curve [l5| (incorpo- 
rated into CORSIKA). This procedure was repeated for 
all showers in the Livni library with the results similar 
to those shown in Figs. Q][3J 

Figure [T] shows the distribution of the reconstructed 
5(600) for showers with thinning simulated with the same 
initial random seed (and thus the same first interaction) 
as three representative e = Livni showers. Though 
quite wide for e = 10~ 4 thinning, the distributions of 
5(600)/5(600) no thinning are well centered at unity. 

The distribution of the mean values of 
5(600)/5(600) no thinning for the ensembles of the 
thinned showers is presented in Fig. [3 for a uniform 
sample of twenty different e = showers. For each 
of them, 500 showers with e = 10~ 4 were simulated 
with the same first interaction as the corresponding 
e = shower. The values of the observable averaged 
over 500 thinned showers approximate the "exact" 
5 (600) no thinning with the accuracy of about 3%, which 
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FIG. 1: Distribution of S(600)/5(600) no thummg, where 
5(600) is reconstructed from 500 showers simulated with 
e = 10" 4 and the same random seed as the corresponding 
e — shower, for three different random seeds (three his- 
tograms). The three showers are all vertical, induced by 
10 18 eV protons at the AGASA location. 
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FIG. 2: Distribution of 5(600 ) average /S (600) no thinning, where 
S'(600)avGragc is the average of the reconstructed 5(600) over 
a sample of 500 showers simulated with e = 10~ 4 and the 
same random seed as the corresponding e — shower, for 
20 different random seeds. The 20 showers are all vertical, 
induced by 10 17 eV protons at the Telescope-Array location. 



is consistent with the level of statistical fluctuations, 
1/V500 ~ 4%. We have found the same distributions 
for other observables considered, p p (1000) and X max . 
The important conclusion is that for the first time, 
the usual assumption that thinning does not introduce 
systematic errors in the reconstructed observables has 
been checked by explicit comparison of e = shower 
and averaged e ^ showers, at least for energies up to 
10 18 eV, observables 5(600), p M (1000) and A max , and 
proton, photon and iron primaries. 

The spread of observables reconstructed from thinned 
showers depends on the thinning level e. It is not the 
width of the distribution but the average deviation of 
the observables from those of an e = shower which is 
the most interesting for practical purposes. This quantity 
is plotted in Fig. [3] for a typical shower from the Livni 
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FIG. 3: RMS deviations from unity of 

5'(600)/S'(600) 11 o thinning (blue triangles and thick blue 
line), p M (1000)/p,j(1000)no thinning (red boxes and red dotted 

line), X ma x/X m ax ! no thinning 

(green diamonds and green 
dashed line), where quantities with subscript "no thinning" 
are reconstructed from a 10 18 eV proton shower with zenith 
angle 45°, simulated without thinning for the AGASA 
observational conditions, while the rest of quantities are 
reconstructed from large samples of showers simulated with 
various thinning levels for the same input parameters and the 
same initial random seed. "—5 opt" denotes e = 10~ 5 with 
weights limitation; "MS" denotes multisampling (20 x 10 -4 ) 
discussed in Sec. IIVI 
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library. 

We note in passing that, technically, to study the 
spread at a given e with CORSIKA, one has to simu- 
late showers with slightly different thinning levels (other- 
wise they all would be absolutely identical, given a fixed 
random seed). For instance, to obtain the points cor- 
responding to e = 10~ 5 in Fig. [31 we have simulated 
500 showers with different thinning levels in the interval 
0.99 • 10~ 5 < e < 1.01 • 10~ 5 . 



B. Distributions of showers 

In most cases one is not interested in details of a par- 
ticular realization of a shower; it is the ensemble of sim- 
ulated showers with fixed initial parameters but varied 
random seeds which is compared to the real data. The 
study of Sec. IIII Al does not help seriously to estimate 
the effect of thinning on these distributions of parame- 
ters because the size of fluctuations seen, e.g., in Fig. [3] 
is determined by a combination of artificial fluctuations 
and a part of real ones: the random seed together with 
initial conditions fixes the first interaction, but different 
thinning levels introduce variations in other interactions 
and effectively change the simulation of the entire shower 
development. 

To estimate the effect of thinning on the distribution of 
observables, we have simulated samples of showers with 
fixed initial conditions but different random seeds for var- 
ious thinning levels, including e = 0. We have considered 
samples of E = 10 17 eV vertical proton-induced showers 



FIG. 4: Width of the S(600) distribution (upper panel) and 
p M (1000) distribution (lower panel) for 10 17 eV vertical proton 
showers simulated with and without thinning for the Tele- 
scope Array observational conditions. Statistical error bars 
are due to a limited number of simulated showers. 



consisting of 20 showers with e = 0, 100 showers with 
e = 10~ 5 , 100 showers with e = 10 -5 and weight limita- 
tion, 100 showers with e = 10~ 4 and 100 showers with 
e = 10~ 4 and weight limitation. 

The simulations have been performed using 
QGSJET II and GHEISHA as hadronic interaction 
models, for the observational conditions of the Telescope 
Array experiment. The distributions of S^OO), p M (1000) 
and X max have been reconstructed with statistical 
fluctuations (originated from the limited number n 
of showers in the samples) of about 1/^/n, that is, 
about 23% for e = showers and about 10% for the 
other samples. Figure [4] illustrates the widths of the 
distributions obtained at different e. Artificial fluctu- 
ations in 5(600) and p M (1000) caused by thinning are 
clearly seen by comparing e = 10~ 4 case with others 
(for X max the artificial fluctuations are quite small). 
We note that, for a given e, the fluctuations should be 
larger at high energy since the multiplicity of hadronic 
interactions grows with energy and thinning starts to 
operate earlier in the shower affecting the first few 
interactions which determine the fluctuations. For many 
practical purposes, these artificial fluctuations should be 
efficiently suppressed. 
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IV. MULTISAMPLING: AN ECONOMICAL 
METHOD TO SUPPRESS ARTIFICIAL 
FLUCTUATIONS 

From the results of the previous section, we conclude 
that the use of thinning is well motivated when one is 
interested in the reconstruction of the central values of 
fluctuating observables (the most important application 
is e.g. to establish a relation between, say, 5(600) and en- 
ergy for a given experimental setup). On the other hand, 
thinning may limit the precision of composition studies, 
where the observed value of some quantity is compared 
to the simulated distributions of the same quantity for 
different primaries, and the width of these distributions 
is of crucial importance (see e.g. the proton-iron com- 
parison in examples of Ref. [161]). 

As it has been pointed out above, the effect of physical 
fluctuations on the distribution of an observable quantity 
should be, in principle, estimated by simulating a set of 
showers with the same physical parameters, with differ- 
ent random seeds and without thinning. To obtain a good 
approximation to this distribution, we make use of the 
results of Sec. IIII Al (see, in particular, Figs. [1] and 
The average of an observable over a sample of thinned 
showers with a fixed initial random seed approximates 
the value of the same observable for an e = shower 
with the same random seed with a good accuracy. The 
distribution of observables for e = showers with differ- 
ent random seeds is then approximated by a distribution 
of these approximated observables calculated for samples 
with random seeds varying from one sample to another 
but fixed inside a sample. A practical way to do this is 
as follows: 

• instead of a single shower with e = 0, simulate N 
showers with some e = eo ^ and fixed random 
seed; 

• reconstruct the observable for each of N showers, 
average over these N realizations and keep this av- 
erage value which approximates the result for a sin- 
gle shower without thinning; 

• repeat the procedure M times for different random 
seeds to mimic a simulation of M showers without 
thinning and to obtain the required distribution of 
the observable. 

We will refer to this procedure as to multisampling 
(N x eo). Even for relatively large e, averaging over a 
sufficiently large number of showers (JV) gives a good ap- 
proximation to an eo = value of the observable; the 
larger N, the better the approximation. The required 
value of N may be estimated as follows. Consider the 
distribution of an observable reconstructed from show- 
ers simulated with the thinning level close to eo f° r a 
given initial random seed. Assume that the distribution 
is Gaussian with the width a (though the qualitative con- 
clusions do not depend on the exact form of the distri- 
bution, we note that, in practice, it is indeed very close 
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FIG. 5: Normalized distributions of the reconstructed 5(600) 
(in VEM/m 2 ): dashed line, 20 showers without thinning; solid 
line, 500 showers with (20 x 10 -4 ) multisampling. All showers 
are vertical, induced by 10 17 eV primary protons simulated 
for the conditions of the Telescope Array. Vertical error bars 
are due to limited statistics. 



to Gaussian [T7|); then one needs N measurements to 
know the mean value with the precision ~ a/yN. Nu- 
merical results for the Livni showers demonstrate that 
(N x eo) multisampling for N ~ 15 ... 20 and eo ~ 10~ 4 
results in the precision of 3 4- 4% in the reconstruction of 
5(600), p M (1000) and A" max of the original e = show- 
ers. The distributions of parameters reconstructed from 
showers without thinning are consistent (within statisti- 
cal errors) with those extracted by making use of multi- 
sampling. The distributions of 5(600) are presented in 
Fig. El 

In Fig. [5] we present the widths of the distributions ob- 
tained with the usual thinning and with multisampling 
for E = 5 • 10 19 eV vertical proton-induced showers; the 
limited statistics (we used n = 200 showers) implies the 
statistical uncertainty of about l/y/n ~ 7%. The gain in 
precision is clearly seen; for the case of 5 TO 19 eV the mul- 
tisampled distribution (which is expected to mimic the 
e = distribution with a good accuracy) allows us to es- 
timate the size of purely artificial fluctuations caused by 
thinning. For instance, for e = 10~ 5 with weights limita- 
tions, these fluctuations remain at the level of > 10% for 
5(600) and of > 12% for ^(1000). Let us note in pass- 
ing that, for this particular simulation (5 • 10 19 eV verti- 
cal protons at the Telescope Array location) and for our 
choice of hadronic models (QGSJET II and GHEISHA), 
the choice of maximal weights suggested in Ref. [ij may 
not be optimal. 

Let us compare now the computer resources needed 
for calculations with the standard thinning (with and 
without weights limitations) and with multisampling. 

The disk space scales as the number of simulated par- 
ticles; Fig. [7] illustrates this fact. We see that the multi- 
sampling (20 x 10~ 4 ) saves the disk space compared to 
e = 10~ 5 with weights limitation, giving at the same time 



6 



S(600) 



10" b 
optimal 10" 5 

10- 4 

MS 20x1 0" 4 
MS 10x10" 4 





a. 


- ±_ 


+ 






+ 







-15% -5% 5% 15% 



p„(1000) 



10" 



optimal 10 



10" 



MS 20x10" 



MS 10x10" 











+ 






-i 




+ 















-30% -15% 



15% 30% 



FIG. 6: Width of the 5(600) distribution (upper panel) and 
p M (1000) distribution (lower panel) for 200 showers initiated 
by 5 ■ 10 19 eV vertical protons simulated with thinning and 
with multisampling for the Telescope Array observational con- 
ditions. Statistical error bars are due to a limited number of 
simulated showers. The choice of maximal weights suggested 
in Ref. [J] may not be optimal for this study (see the text). 
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FIG. 7: The size of the CORSIKA output file for different 
thinning levels and procedures (a 10 18 eV proton shower with 
zenith angle of 45°, AG AS A observational conditions). 



gain in the precision of simulations. 

The CPU time is very sensitive to the choice of the 
hadronic interaction model: since thinning starts to work 
when the number of particles is large enough, the first 
few interactions are simulated in full even for relatively 
large e. If the high-energy model is slow, then the ef- 
fect of multisampling on the computational time is not 
so pronounced. By variations of the hadronic interaction 
models, wc have estimated the average time consumed by 
QGSJET II, SYBILL, FLUKA and GHEISHA for sim- 
ulations of showers at energies 10 17 eV and 5 • 10 19 eV. 
For 5 • 10 19 eV vertical proton showers, (20 x 10 -4 ) mul- 
tisampling is about 5 times faster than 10 -5 thinning 
with weights limitation for SYBILL while for (very slow) 
QGSJET II, both take roughly the same time. A way to 
change the multisampling procedure in order to gain in 
the CPU time for any hadronic model is discussed below 
in Sec. fVl 



V. DISCUSSION AND CONCLUSIONS 

A library of atmospheric showers has been simulated 
without the thinning approximation. The showers have 
been used for a quantitative direct study of the effect of 
thinning on the reconstruction of signal (S) and muon 
(p M ) densities at the ground level as well as on the depth 
A'max of the maximal shower development. We demon- 
strate that thinning does not introduce systematic shifts 
into these observables, as was conjectured but never ex- 
plicitly checked. We estimate the size of artificial fluctu- 
ations which appear due to the reduction of the number 
of particles in the framework of the thinning approxima- 
tion; these unphysical fluctuations may affect the preci- 
sion, e.g., of the composition studies. For instance, at 
the energies of 5 • 10 19 eV for vertical proton primaries, 
artificial fluctuations are about 10% in the signal density 
at 600 m and about 12% in the muon density at 1000 m 
for e = 10~ 5 thinning with weight limitations. An effec- 
tive method to suppress these artificial fluctuations, mul- 
tisampling, is suggested and studied. The method does 
not invoke any changes in simulation codes; only parame- 
ters of, say, the CORSIKA input are affected. Compared 
to the 10~ 5 thinning with weights limitations, it gives 
a similar precision but allows one to gain an order-of- 
magnitude decrease in the required disk space. Gain in 
the CPU time depends on the speed of the high-energy 
interaction model: it is of order 5 10 for fast ones 
(SYBILL) and of order 1 for slow ones (QGSJET II). 

A way to change the multisampling procedure in order 
to further improve the gain in the CPU time is to simulate 
the high-energy part of a shower once for each initial ran- 
dom seed while having the low-energy part multisampled. 
The multisampling procedure described above is a partic- 
ular case of such improved procedure with a high-energy 
part restricted to the first interaction only. We would 
expect the modification to make it possible to conserve 
the physical fluctuations in the second and several follow- 
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ing interactions and will allow for an order-of-magnitude 
improvement in the computational time for any hadronic 
model. However, it would require (simple) changes in the 
simulation codes thus loosing an important advantage 
of the multisampling discussed above: to implement the 
latter, one operates with the standard simulation code 
(e.g., CORSIKA) without any modifications. This mini- 
mal change is to add the option to start simulations from 
a predefined set of the primary particles. 
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